The following function can be used to exchange, or permute, the
elements of a vector.
- Function: int gsl_vector_swap_elements (gsl_vector * V, size_t i,
size_t j)
This function exchanges the I-th and J-th elements of the vector V
in-place.
- Function: int gsl_vector_reverse (gsl_vector * V)
This function reverses the order of the elements of the vector V.
File: gsl-ref.info, Node: Vector operations, Next: Finding maximum and minimum elements of vectors, Prev: Exchanging elements, Up: Vectors and Matrices
Vector operations
=================
- Function: int gsl_vector_add (gsl_vector * A, const gsl_vector * B)
This function adds the elements of vector B to the elements of
vector A, a'_i = a_i + b_i. The two vectors must have the same
length.
- Function: int gsl_vector_sub (gsl_vector * A, const gsl_vector * B)
This function subtracts the elements of vector B from the elements
of vector A, a'_i = a_i - b_i. The two vectors must have the same
length.
- Function: int gsl_vector_mul (gsl_vector * A, const gsl_vector * B)
This function multiplies the elements of vector A by the elements
of vector B, a'_i = a_i * b_i. The two vectors must have the same
length.
- Function: int gsl_vector_div (gsl_vector * A, const gsl_vector * B)
This function divides the elements of vector A by the elements of
vector B, a'_i = a_i / b_i. The two vectors must have the same
length.
- Function: int gsl_vector_scale (gsl_vector * A, const double X)
This function multiplies the elements of vector A by the constant
factor X, a'_i = x a_i.
- Function: int gsl_vector_add_constant (gsl_vector * A, const double
X)
This function adds the constant value X to the elements of the
vector A, a'_i = a_i + x.
File: gsl-ref.info, Node: Finding maximum and minimum elements of vectors, Next: Vector properties, Prev: Vector operations, Up: Vectors and Matrices
This function returns the indices of the minimum and maximum
values in the vector V, storing them in IMIN and IMAX. When there
are several equal minimum or maximum elements then the lowest
indices are returned.
File: gsl-ref.info, Node: Vector properties, Next: Example programs for vectors, Prev: Finding maximum and minimum elements of vectors, Up: Vectors and Matrices
Vector properties
=================
- Function: int gsl_vector_isnull (const gsl_vector * V)
This function returns 1 if all the elements of the vector V are
zero, and 0 otherwise.
File: gsl-ref.info, Node: Example programs for vectors, Next: The matrix struct, Prev: Vector properties, Up: Vectors and Matrices
Example programs for vectors
============================
This program shows how to allocate, initialize and read from a vector
using the functions `gsl_vector_alloc', `gsl_vector_set' and
`gsl_vector_get'.
#include <stdio.h>
#include <gsl/gsl_vector.h>
int main ()
{
int i;
gsl_vector * v = gsl_vector_alloc (3) ;
for (i = 0; i < 3; i++)
{
gsl_vector_set (v, i, 1.23 + i);
}
for (i = 0; i < 100; i++)
{
printf("v_%d = %g\n", i, gsl_vector_get (v, i));
}
}
Here is the output from the program. The final loop attempts to read
outside the range of the vector `v', and the error is trapped by the
range-checking code in `gsl_vector_get'.
v_0 = 1.23
v_1 = 2.23
v_2 = 3.23
gsl: vector_source.c:12: ERROR: index out of range
IOT trap/Abort (core dumped)
The next program shows how to write a vector to a file.
#include <stdio.h>
#include <gsl/gsl_vector.h>
int main ()
{
int i;
gsl_vector * v = gsl_vector_alloc (100) ;
for (i = 0; i < 100; i++)
{
gsl_vector_set (v, i, 1.23 + i);
}
{
FILE * f = fopen("test.dat", "w") ;
gsl_vector_fprintf (f, v, "%.5g");
fclose (f);
}
}
After running this program the file `test.dat' should contain the
elements of `v', written using the format specifier `%.5g'. The vector
could then be read back in using the function `gsl_vector_fscanf (f,
v)'.
File: gsl-ref.info, Node: The matrix struct, Next: Matrix allocation, Prev: Example programs for vectors, Up: Vectors and Matrices
The matrix struct
=================
Matrices are defined by a `gsl_matrix' structure which describes a
generalized slice of a block. Like a vector it represents a "view" of
an area of memory, but uses two indices instead of one.
The `gsl_matrix' structure contains five members, the two dimensions
of the matrix, a physical dimension, a pointer to the memory where the
elements of the matrix are stored, DATA, and a pointer to the block
owned by the matrix BLOCK, if any. The physical dimension determines
the memory layout and can differ from the matrix dimension to allow the
use of submatrices. The `gsl_matrix' structure is very simple and
looks like this,
typedef struct
{
size_t size1;
size_t size2;
size_t tda;
double * data;
gsl_block * block;
} gsl_matrix ;
Matrices are stored in row-major order, meaning that each row of
elements forms a contiguous block in memory. This is the standard
"C-ordering" of two-dimensional arrays. Note that FORTRAN stores arrays
in column-major order. The number of rows is SIZE1. The range of valid
row indices runs from 0 to `size1-1'. Similarly SIZE2 is the number of
columns. The range of valid column indices runs from 0 to `size2-1'.
The physical row dimension TDA, or "trailing dimension", specifies the
size of a row of the matrix as laid out in memory.
For example, in the following matrix SIZE1 is 3, SIZE2 is 4, and TDA
is 8. The physical memory layout of the matrix begins in the top left
hand-corner and proceeds from left to right along each row in turn.
00 01 02 03 XX XX XX XX
10 11 12 13 XX XX XX XX
20 21 22 23 XX XX XX XX
Each unused memory location is represented by "`XX'". The pointer DATA
gives the location of the first element of the matrix in memory. The
pointer BLOCK stores the location of the memory block owned by the
matrix (if any). This block will be deallocated when the matrix is
freed. If the matrix is only a view of another object then the pointer
BLOCK is null.
The functions for allocating and accessing matrices are defined in
`gsl_matrix.h'
File: gsl-ref.info, Node: Matrix allocation, Next: Accessing matrix elements, Prev: The matrix struct, Up: Vectors and Matrices
Matrix allocation
=================
The functions for allocating memory to a matrix follow the style of
`malloc' and `free'. They also perform their own error checking. If
there is insufficient memory available to allocate a vector then the
functions call the GSL error handler (with an error number of
`GSL_ENOMEM') in addition to returning a null pointer. Thus if you use
the library error handler to abort your program then it isn't necessary
where the index I runs from 0 to `n1-1' and the index J runs from
0 to `n2-1'.
The `data' pointer of the returned matrix struct is set to null if
the combined parameters (K1,K2,N1,N2,TDA) overrun the ends of the
original matrix.
The new matrix is only a view of the block underlying the existing
matrix, M. The block containing the elements of M is not owned by
the new matrix. When the new matrix goes out of scope the
original matrix M and its block will continue to exist. The
original memory can only be deallocated by freeing the original
matrix. Of course, the original matrix should not be deallocated
while the new matrix is still in use.
File: gsl-ref.info, Node: Creating row and column views, Next: Copying matrices, Prev: Creating submatrix views, Up: Vectors and Matrices
Creating row and column views
=============================
In general there are two ways to access an object, by reference or by
copying. The functions described in this section create vectors which
allow access to a row or column of a matrix by reference. Modifying
elements of the vector is equivalent to modifying the matrix, since both
the vector and the matrix point to the same memory block.
- Function: gsl_vector gsl_matrix_row (gsl_matrix * M, size_t I)
This function returns `gsl_vector' struct which is a view of the
I-th row of the matrix M. The `data' pointer of the new vector is
set to null if I is out of range.
- Function: gsl_vector gsl_matrix_column (gsl_matrix * M, size_t J)
This function returns a `gsl_vector' struct which is a view of the
J-th column of the matrix M. The `data' pointer of the new vector
is set to null if J is out of range.
- Function: gsl_vector gsl_matrix_diagonal (gsl_matrix * M)
This function returns a `gsl_vector' struct which is a view of the
diagonal of the matrix M. The matrix M is not required to be
square. For a rectangular matrix the length of the diagonal is the
same as the smaller dimension of the matrix.
File: gsl-ref.info, Node: Copying matrices, Next: Copying rows and columns, Prev: Creating row and column views, Up: Vectors and Matrices
Copying matrices
================
- Function: int gsl_matrix_memcpy (gsl_matrix * DEST, const gsl_matrix
* SRC)
This function copies the elements of the matrix SRC into the
matrix DEST. The two matrices must have the same size.
- Function: int gsl_matrix_swap (gsl_matrix * M1, const gsl_matrix *
M2)
This function exchanges the elements of the matrices M1 and M2 by
copying. The two matrices must have the same size.
File: gsl-ref.info, Node: Copying rows and columns, Next: Exchanging rows and columns, Prev: Copying matrices, Up: Vectors and Matrices
Copying rows and columns
========================
The functions described in this section copy a row or column of a
matrix into a vector. This allows the elements of the vector and the
matrix to be modified independently. Note that if the matrix and the
vector point to overlapping regions of memory then the result will be
undefined.
- Function: int gsl_matrix_get_row (gsl_vector * V, const gsl_matrix *
M, size_t I)
This function copies the elements of the I-th row of the matrix M
into the vector V. The length of the vector must be the same as
the length of the row.
- Function: int gsl_matrix_get_col (gsl_vector * V, const gsl_matrix *
M, size_t J)
This function copies the elements of the I-th column of the matrix
M into the vector V. The length of the vector must be the same as
the length of the column.
- Function: int gsl_matrix_set_row (gsl_matrix * M, size_t I, const
gsl_vector * V)
This function copies the elements of the vector V into the I-th
row of the matrix M. The length of the vector must be the same as
the length of the row.
- Function: int gsl_matrix_set_col (gsl_matrix * M, size_t J, const
gsl_vector * V)
This function copies the elements of the vector V into the I-th
column of the matrix M. The length of the vector must be the same
as the length of the column.
File: gsl-ref.info, Node: Exchanging rows and columns, Next: Matrix operations, Prev: Copying rows and columns, Up: Vectors and Matrices
Exchanging rows and columns
===========================
The following functions can be used to exchange the rows and columns
of a matrix.
- Function: int gsl_matrix_swap_rows (gsl_matrix * M, size_t I, size_t
J)
This function exchanges the I-th and J-th rows of the matrix M
in-place.
- Function: int gsl_matrix_swap_columns (gsl_matrix * M, size_t I,
size_t J)
This function exchanges the I-th and J-th columns of the matrix M
in-place.
- Function: int gsl_matrix_swap_rowcol (gsl_matrix * M, size_t I,
size_t J)
This function exchanges the I-th row and J-th column of the matrix
M in-place. The matrix must be square for this operation to be
possible.
- Function: int gsl_matrix_transpose (gsl_matrix * M)
This function replaces the matrix M by its transpose by copying
the elements of the matrix in-place. The matrix must be square
for this operation to be possible.
File: gsl-ref.info, Node: Matrix operations, Next: Finding maximum and minimum elements of matrices, Prev: Exchanging rows and columns, Up: Vectors and Matrices
Matrix operations
=================
- Function: int gsl_matrix_add (gsl_matrix * A, const gsl_matrix * B)
This function adds the elements of matrix B to the elements of
matrix A, a'(i,j) = a(i,j) + b(i,j). The two matrices must have the
same dimensions.
- Function: int gsl_matrix_sub (gsl_matrix * A, const gsl_matrix * B)
This function subtracts the elements of matrix B from the elements
of matrix A, a'(i,j) = a(i,j) - b(i,j). The two matrices must have
the same dimensions.
- Function: int gsl_matrix_mul (gsl_matrix * A, const gsl_matrix * B)
This function multiplies the elements of matrix A by the elements
of matrix B, a'(i,j) = a(i,j) * b(i,j). The two matrices must have
the same dimensions.
- Function: int gsl_matrix_div (gsl_matrix * A, const gsl_matrix * B)
This function divides the elements of matrix A by the elements of
matrix B, a'(i,j) = a(i,j) / b(i,j). The two matrices must have the
same dimensions.
- Function: int gsl_matrix_scale (gsl_matrix * A, const double X)
This function multiplies the elements of matrix A by the constant
factor X, a'(i,j) = x a(i,j).
- Function: int gsl_matrix_add_constant (gsl_matrix * A, const double
X)
This function adds the constant value X to the elements of the
matrix A, a'(i,j) = a(i,j) + x.
File: gsl-ref.info, Node: Finding maximum and minimum elements of matrices, Next: Matrix properties, Prev: Matrix operations, Up: Vectors and Matrices
Finding maximum and minimum elements of matrices
================================================
- Function: double gsl_matrix_max (const gsl_matrix * M)
This function returns the maximum value in the matrix M.
- Function: double gsl_matrix_min (const gsl_matrix * M)
This function returns the minimum value in the matrix M.
- Function: void gsl_matrix_minmax (const gsl_matrix * M, double *
MIN_OUT, double * MAX_OUT)
This function returns the minimum and maximum values in the matrix
M, storing them in MIN_OUT and MAX_OUT.
- Function: void gsl_matrix_max_index (const gsl_matrix * M, size_t *
IMAX, size_t * JMAX)
This function returns the indices of the maximum value in the
matrix M, storing them in IMAX and JMAX. When there are several
equal maximum elements then the first element found is returned.
- Function: void gsl_matrix_min_index (const gsl_matrix * M, size_t *
IMAX, size_t * JMAX)
This function returns the indices of the minimum value in the
matrix M, storing them in IMAX and JMAX. When there are several
equal minimum elements then the first element found is returned.
- Function: void gsl_matrix_minmax_index (const gsl_matrix * M, size_t
* IMIN, size_t * IMAX)
This function returns the indices of the minimum and maximum
values in the matrix M, storing them in (IMIN,JMIN) and
(IMAX,JMAX. When there are several equal minimum or maximum
elements then the first elements found are returned.
File: gsl-ref.info, Node: Matrix properties, Next: Example programs for matrices, Prev: Finding maximum and minimum elements of matrices, Up: Vectors and Matrices
Matrix properties
=================
- Function: int gsl_matrix_isnull (const gsl_matrix * M)
This function returns 1 if all the elements of the matrix M are
zero, and 0 otherwise.
File: gsl-ref.info, Node: Example programs for matrices, Next: Vector and Matrix References and Further Reading, Prev: Matrix properties, Up: Vectors and Matrices
Example programs for matrices
=============================
This program shows how to allocate, initialize and read from a matrix
using the functions `gsl_matrix_alloc', `gsl_matrix_set' and
`gsl_matrix_get'.
#include <stdio.h>
#include <gsl/gsl_matrix.h>
int main ()
{
int i, j;
gsl_matrix * m = gsl_matrix_alloc (10, 3) ;
for (i = 0; i < 10; i++)
for (j = 0; j < 3; j++)
gsl_matrix_set (m, i, j, 0.23 + 100*i + j);
for (i = 0; i < 100; i++)
for (j = 0; j < 3; j++)
printf("m(%d,%d) = %g\n", i, j, gsl_matrix_get (m, i, j));
}
Here is the output from the program. The final loop attempts to read
outside the range of the matrix `m', and the error is trapped by the
range-checking code in `gsl_matrix_get'.
m(0,0) = 0.23
m(0,1) = 1.23
m(0,2) = 2.23
m(1,0) = 100.23
m(1,1) = 101.23
m(1,2) = 102.23
...
m(9,2) = 902.23
gsl: matrix_source.c:13: ERROR: first index out of range
IOT trap/Abort (core dumped)
The next program shows how to write a matrix to a file.
#include <stdio.h>
#include <gsl/gsl_matrix.h>
int main ()
{
int i, j, differences = 0;
gsl_matrix * m = gsl_matrix_alloc (100, 100) ;
gsl_matrix * a = gsl_matrix_alloc (100, 100) ;
for (i = 0; i < 100; i++)
for (j = 0 ; j < 100; j++)
gsl_matrix_set (m, i, j, 0.23 + i + j);
{
FILE * f = fopen("test.dat", "w") ;
gsl_matrix_fwrite (f, m);
fclose (f);
}
{
FILE * f = fopen("test.dat", "r") ;
gsl_matrix_fread (f, a);
fclose (f);
}
for (i = 0; i < 100; i++)
for (j = 0 ; j < 100; j++)
if (gsl_matrix_get(m, i, j) != gsl_matrix_get(a, i, j))
differences ++ ;
printf("differences = %d (should be zero)\n", differences) ;
}
After running this program the file `test.dat' should contain the
elements of `m', written in binary format. The matrix which is read
back in using the function `gsl_matrix_fread' should be exactly equal
to the original matrix.
File: gsl-ref.info, Node: Vector and Matrix References and Further Reading, Prev: Example programs for matrices, Up: Vectors and Matrices
References and Further Reading
==============================
The block, vector and matrix objects in GSL follow the `valarray'
model of C++. A description of this model can be found in the following
reference,
B. Stroustrup, `The C++ Programming Language' (3rd Ed), Section
22.4 Vector Arithmetic. Addison-Wesley 1997, ISBN 0-201-88954-4.
File: gsl-ref.info, Node: BLAS Support, Next: Linear Algebra, Prev: Vectors and Matrices, Up: Top
BLAS Support
************
* Menu:
* Introduction::
* Organization::
* GSL BLAS Interface::
* Raw BLAS Interface::
* BLAS References and Further Reading::
File: gsl-ref.info, Node: Introduction, Next: Organization, Up: BLAS Support
Introduction
============
The Basic Linear Algebra Standard (BLAS) prescribes a set of
fundamental operations on vectors and matrices which can be used to
create higher-level linear algebra functionality. Because of the many
possible types of vectors and matrices, including variations in base
precision, the number of functions specified by the BLAS standard is
tediously large.
Furthermore, because of the long history of BLAS and BLAS-like
implementations, and the need for interoperability, the complexity of
this part of GSL is somewhat high. Casual users will often find what
they need in the `gsl_linalg' library, which implements a set of common
high-level linear algebra operations on GSL vector and matrix objects.
Some familiarity with the BLAS standards, including both the legacy and
draft interface standards, is assumed in this chapter.
